In [1]:
%matplotlib qt
In [2]:
import os.path as op
import numpy as np
import matplotlib.pyplot as plt
import operator
import mne
In [3]:
data_path = "/Users/nielskjaergaardmadsen/Documents/CFIN/Data/bmt_cleaned/run_3_raw_sss.fif"
raw = mne.io.read_raw_fif(data_path)
In [4]:
events = mne.find_events(raw, min_duration=0.01)
In [5]:
event_id = {'Start':10,
'Q1/1':101,'Q1/2':102,'Q1/3':103,'Q2/1':104,'Q2/2':105,'Q2/3':106,
'Q3/1':107,'Q3/2':108,'Q3/3':109,'Q4/1':110,'Q4/2':111,'Q4/3':112,
'Q1/4':113,'Q1/5':114,'Q1/6':115,'Q2/4':116,'Q2/5':117,'Q2/6':118,
'Q3/4':119,'Q3/5':120,'Q3/6':121,'Q4/4':122,'Q4/5':123,'Q4/6':124,
'Q1/7':125,'Q1/8':126,'Q1/9':127,'Q2/7':128,'Q2/8':129,'Q2/9':130,
'Q3/7':131,'Q3/8':132,'Q3/9':133,'Q4/7':134,'Q4/8':135,'Q4/9':136}
In [6]:
l_freq, h_freq = 1,40
raw.load_data()
raw.filter(l_freq, h_freq)
Out[6]:
In [ ]:
raw.plot(order='selection')
In [ ]:
raw.plot?
In [7]:
tmin, tmax = -0.05, 0.2
In [14]:
baseline = None
reject = {'mag':4e-12}
#epochs = mne.Epochs(raw, events=events,baseline = baseline, event_id=event_id, tmin=tmin,
# tmax=tmax, reject=reject, preload=True,detrend = 0)
In [ ]:
epochs.plot(block=True)
In [ ]:
epochs.plot_drop_log()
In [ ]:
# Q1
ts_args = dict(gfp=True)
topomap_args = dict(sensors=False)
x = 0
fig, ax = plt.subplots(9,4)
event_Q1_id = {'Q1/1':101,'Q1/2':102,'Q1/3':103,'Q1/4':113,'Q1/5':114,'Q1/6':115,
'Q1/7':124,'Q1/8':125,'Q1/9':126}
times = [0.035, 0.085, 0.135, 0.2] # 0.2 is not in use
for key in event_Q1_id.keys():
evoked = epochs[key].average()
#if x==8:
print(key)
evoked.plot_topomap(times = times,axes=ax[x], show=True)
#evoked.plot(exclude=(),axes = ax2[x])
#evoked.plot_joint(title=key, times=[0.035, 0.085, 0.135],ts_args=ts_args, topomap_args=topomap_args)
#else:
#evoked.plot_topomap(times = times[0],axes=ax[x], show=False)
#evoked.plot_topomap(times = times[1],axes=ax2[x], show=False)
#evoked.plot(exclude=(),axes = ax2[x])
x +=1
print('Calculating')
if x==9:
print('Breaking')
break
#evoked_Q1_4 = epochs['Q1/2'].average()
#evoked_Q1_7 = epochs['Q1/3'].average()
In [ ]:
# Q2
ts_args = dict(gfp=True)
topomap_args = dict(sensors=False)
x = 0
fig, ax = plt.subplots(9,4)
event_Q2_id = {'Q2/1':104,'Q2/2':105,'Q2/3':106,
'Q2/4':116,'Q2/5':117,'Q2/6':118,
'Q2/7':127,'Q2/8':128,'Q2/9':129}
times = [0.035, 0.085, 0.135, 0.2] # 0.2 is not in use
for key in event_Q2_id.keys():
evoked = epochs[key].average()
#if x==8:
evoked.plot_topomap(times = times,axes=ax[x], show=True)
#evoked.plot(exclude=(),axes = ax2[x])
#evoked.plot_joint(title=key, times=[0.035, 0.085, 0.135],ts_args=ts_args, topomap_args=topomap_args)
#else:
#evoked.plot_topomap(times = times[0],axes=ax[x], show=False)
#evoked.plot_topomap(times = times[1],axes=ax2[x], show=False)
#evoked.plot(exclude=(),axes = ax2[x])
x +=1
print('Calculating')
if x==9:
print('Breaking')
break
#evoked_Q1_4 = epochs['Q1/2'].average()
#evoked_Q1_7 = epochs['Q1/3'].average()
In [ ]:
for var in sorted_:
a = var
print(a)
In [ ]:
evoked.plot?
In [ ]:
ts_args = dict(gfp=True)
topomap_args = dict(sensors=False)
sorted_ = sorted(event_id.items(), key=operator.itemgetter(1))
sorted_.remove(('Start', 10))
times = [0.085]
fig, ax = plt.subplots(9,1)
fig2, ax2 = plt.subplots(3,2)
plt.tight_layout()
x = 0
y = 0
quadrant = 0 # 0 for Q1, 3 for Q2, 6 for Q3 and 9 for Q4
for i in range(0, 3):
for ii in range(0, 3):
if i > 0:
ind = ii+6*y+quadrant
print(ind)
key = sorted_[ind]
print(key[0])
evoked = epochs[key[0]].average()
#evoked.plot(gfp=True,axes = ax2[x])
evoked.plot_topomap(times = times,axes=ax[x], show=False,colorbar = False, title='Q1')
x +=1
else:
ind = ii+quadrant
print(ind)
key = sorted_[ind]
print(key[0])
evoked = epochs[key[0]].average()
evoked.plot(gfp=True,axes = ax2[x])
evoked.plot_topomap(times = times,axes=ax[x], show=False,colorbar = False)
x +=1
y += 2
In [ ]:
times = np.arange(0.05, 0.15, 0.01)
evoked.plot_topomap(times, ch_type='grad')
fig, anim = evoked.animate_topomap(ch_type='mag', times=times, frame_rate=10)
evoked.plot_topomap(0.1, ch_type='grad', show_names=True, colorbar=False,
size=6, res=128, title='Visual respons')
In [ ]:
fig = evoked_Q1_1.plot(exclude=())
fig2 = evoked_Q2_1.plot(exclude=())
fig3 = evoked_Q3_1.plot(exclude=())
fig4 = evoked_Q4_1.plot(exclude=())
In [ ]:
ts_args = dict(gfp=True)
topomap_args = dict(sensors=False)
#evoked_Q1_1.plot_joint(title='Q1/2', times=[0.035, 0.085, 0.135],
#ts_args=ts_args, topomap_args=topomap_args)
#evoked_Q2_1.plot_joint(title='Q2/2', times=[0.035, 0.085, 0.135],
# ts_args=ts_args, topomap_args=topomap_args)
#evoked_Q3_1.plot_joint(title='Q3/2', times=[0.035, 0.085, 0.135],
# ts_args=ts_args, topomap_args=topomap_args)
#evoked_Q4_1.plot_joint(title='Q4/2', times=[0.035, 0.085, 0.135],
# ts_args=ts_args, topomap_args=topomap_args)
evoked_Q1_9 = epochs['Q1/9'].average()
evoked_Q1_9.plot_joint(title='Q1/9', times=[0.035, 0.085, 0.135],
ts_args=ts_args, topomap_args=topomap_args)
In [ ]:
evoked_Q1_1.plot_topomap(times = 'peaks', title = 'Q1/1')
evoked_Q1_4.plot_topomap(times = 'peaks', title = 'Q1/4')
evoked_Q1_7.plot_topomap(times = 'peaks', title = 'Q1/7')
In [ ]:
evoked_Q2_1.plot_topomap(times = 'peaks', title = 'Q2/1')
evoked_Q2_4.plot_topomap(times = 'peaks', title = 'Q2/4')
evoked_Q2_7.plot_topomap(times = 'peaks', title = 'Q2/7')
In [ ]:
evoked_Q3_1.plot_topomap(times = 'peaks', title = 'Q3/1')
evoked_Q4_1.plot_topomap(times = 'peaks', title = 'Q4/1')
In [ ]:
raw.plot_psd(tmax=np.inf, fmax=100)
In [ ]:
epochs.plot_psd(fmin=1., fmax=40.)
In [ ]:
epochs.plot_psd_topomap(ch_type='grad', normalize=True)
In [ ]:
subjects_dir = data_path + '/subjects'
trans_fname = data_path + '/MEG/sample/sample_audvis_raw-trans.fif'
maps = mne.make_field_map(evoked_Q1_1, trans=trans_fname, subject='sample',
subjects_dir=subjects_dir, n_jobs=1)
# explore several points in time
field_map = evoked_Q1_1.plot_field(maps, time=.1)
In [ ]:
raw.plot?
In [9]:
from mne import merge_events
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
import mne
from mne.datasets import sample
from mne.decoding import (SlidingEstimator, GeneralizingEstimator,
cross_val_multiscore, LinearModel, get_coef)
In [74]:
tmp_id = np.arange(101,137,1).reshape(12, 3)
upper_id = tmp_id[[0,3,4,7,8,11],:].ravel()
lower_id= tmp_id[[1,2,5,6,9,10],:].ravel()
upper_id
lower_id
Out[74]:
In [75]:
new_events = events.copy()
y = 0
for x in events:
val = x[2]
if val in upper_id:
new_events[y][2] = 1
elif val in lower_id:
new_events[y][2] = 2
y += 1
In [76]:
tst = mne.merge_events(events, lower_id, 2)
In [77]:
tst[:20]
Out[77]:
In [78]:
new_events
Out[78]:
In [79]:
new_event_id = dict(upper=1, lower=2)
epochs_MVP = mne.Epochs(raw = raw, events = new_events,event_id = new_event_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4)
In [80]:
X = epochs_MVP.get_data() # MEG signals: n_epochs, n_channels, n_times
y = epochs_MVP.events[:, 2] # target: Upper or lower
clf = make_pipeline(StandardScaler(), LogisticRegression())
time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc')
scores = cross_val_multiscore(time_decod, X, y, cv=5, n_jobs=1)
# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)
In [81]:
scores.shape
Out[81]:
In [ ]:
# Plot
fig, ax = plt.subplots()
ax.plot(epochs_MVP.times, scores, label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC') # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding')
plt.show()
# You can retrieve the spatial filters and spatial patterns if you explicitly
# use a LinearModel
clf = make_pipeline(StandardScaler(), LinearModel(LogisticRegression()))
time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc')
time_decod.fit(X, y)
coef = get_coef(time_decod, 'patterns_', inverse_transform=True)
evoked = mne.EvokedArray(coef, epochs_MVP.info, tmin=epochs_MVP.times[0])
evoked.plot_joint(times=np.arange(0.,0.2,.50), title='patterns')
In [ ]:
# define the Temporal Generalization object
time_gen = GeneralizingEstimator(clf, n_jobs=1, scoring='roc_auc')
scores = cross_val_multiscore(time_gen, X, y, cv=5, n_jobs=1)
# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
fig, ax = plt.subplots()
ax.plot(epochs_MVP.times, np.diag(scores), label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC')
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Decoding MEG sensors over time')
plt.show()
# Plot the full matrix
fig, ax = plt.subplots(1, 1)
im = ax.imshow(scores, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization Upper/Lower')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
In [10]:
tmp_id_Q = np.arange(101,137,1).reshape(12, 3)
Q1_id = tmp_id_Q[[0,4,8],:].ravel()
Q2_id = tmp_id_Q[[1,5,9],:].ravel()
Q3_id = tmp_id_Q[[2,6,10],:].ravel()
Q4_id = tmp_id_Q[[3,7,11],:].ravel()
In [11]:
tst = mne.merge_events(events, lower_id, 2)
In [12]:
new_events_Q = mne.merge_events(events, Q1_id, 1)
new_events_Q = mne.merge_events(new_events_Q, Q2_id, 2)
new_events_Q = mne.merge_events(new_events_Q, Q3_id, 3)
new_events_Q = mne.merge_events(new_events_Q, Q4_id, 4)
In [17]:
new_event_Q_id = dict(Q1=1, Q2=2, Q3=3, Q4=4)
epochs_MVP_Q = mne.Epochs(raw = raw, events = new_events_Q,event_id = new_event_Q_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4)
In [ ]:
make_pipeline?
In [18]:
X = epochs_MVP_Q.get_data() # MEG signals: n_epochs, n_channels, n_times
y = epochs_MVP_Q.events[:, 2] # target: Q1, Q2, Q3, Q4
In [10]:
scoring_function1 = 'roc_auc'
scoring_function2 = 'accuracy'
In [19]:
clf = make_pipeline(StandardScaler(), LogisticRegression())
print('CLF done')
# time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc')
time_decod = SlidingEstimator(clf, n_jobs=1, scoring=scoring_function2)
print('Time decode')
scores = cross_val_multiscore(time_decod, X, y, cv=5, n_jobs=2)
print('Mean scores')
# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)
print('Plotting')
# Plot
fig, ax = plt.subplots()
ax.plot(epochs_MVP_Q.times, scores, label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC') # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding')
plt.show()
In [21]:
# You can retrieve the spatial filters and spatial patterns if you explicitly
# use a LinearModel
clf = make_pipeline(StandardScaler(), LinearModel(LogisticRegression()))
time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc')
time_decod.fit(X, y)
coef = get_coef(time_decod, 'patterns_', inverse_transform=True)
evoked = mne.EvokedArray(coef, epochs_MVP_Q.info, tmin=epochs_MVP_Q.times[0])
evoked.plot_joint(times=np.arange(0.,0.2,.50), title='patterns')
In [23]:
# define the Temporal Generalization object
time_gen = GeneralizingEstimator(clf, n_jobs=1, scoring=scoring_function)
scores = cross_val_multiscore(time_gen, X, y, cv=5, n_jobs=1)
# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
fig, ax = plt.subplots()
ax.plot(epochs_MVP_Q.times, np.diag(scores), label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('Accuracy')
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Decoding MEG sensors over time in quadrants')
plt.show()
# Plot the full matrix
fig, ax = plt.subplots(1, 1)
im = ax.imshow(scores, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_Q.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
In [66]:
tmp_id_Q = np.arange(101,137,1).reshape(12, 3)
Q1_id = tmp_id_Q[[0,4,8],:].ravel()
Q2_id = tmp_id_Q[[1,5,9],:].ravel()
In [67]:
new_events_Q12 = mne.merge_events(events, Q1_id, 1)
new_events_Q12 = mne.merge_events(new_events_Q12, Q2_id, 2)
In [68]:
new_event_Q12_id = dict(Q1=1, Q2=2)
epochs_MVP_Q12 = mne.Epochs(raw = raw, events = new_events_Q12,event_id = new_event_Q12_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4)
In [69]:
X = epochs_MVP_Q12.get_data() # MEG signals: n_epochs, n_channels, n_times
y = epochs_MVP_Q12.events[:, 2] # target: Q1, Q2
In [70]:
# define the Temporal Generalization object
time_gen = GeneralizingEstimator(clf, n_jobs=1, scoring=scoring_function)
scores = cross_val_multiscore(time_gen, X, y, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
fig, ax = plt.subplots()
ax.plot(epochs_MVP_Q12.times, np.diag(scores), label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('Accuracy')
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Decoding MEG sensors over time in quadrants')
plt.show()
# Plot the full matrix
fig, ax = plt.subplots(1, 1)
im = ax[x].imshow(scores, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_Q.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
In [82]:
tmp_id_Q = np.arange(101,137,1).reshape(12, 3)
R1_id = tmp_id_Q[[0,1,2,3],:].ravel()
R2_id = tmp_id_Q[[4,5,6,7],:].ravel()
R3_id = tmp_id_Q[[8,9,10,11],:].ravel()
In [83]:
new_events_R = mne.merge_events(events, R1_id, 1)
new_events_R = mne.merge_events(new_events_R, R2_id, 2)
new_events_R = mne.merge_events(new_events_R, R3_id, 3)
In [85]:
new_event_R_id13 = dict(R1=1, R3=3)
epochs_MVP_R13 = mne.Epochs(raw = raw, events = new_events_R,event_id = new_event_R_id13, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4)
new_event_R_id12 = dict(R2=1, R3=2)
epochs_MVP_R12 = mne.Epochs(raw = raw, events = new_events_R,event_id = new_event_R_id12, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4)
new_event_R_id23 = dict(R2=1, R3=3)
epochs_MVP_R23 = mne.Epochs(raw = raw, events = new_events_R,event_id = new_event_R_id23, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4)
In [86]:
X13 = epochs_MVP_R13.get_data() # MEG signals: n_epochs, n_channels, n_times
y13 = epochs_MVP_R13.events[:, 2] # target
X12 = epochs_MVP_R12.get_data() # MEG signals: n_epochs, n_channels, n_times
y12 = epochs_MVP_R12.events[:, 2] # target
X23 = epochs_MVP_R23.get_data() # MEG signals: n_epochs, n_channels, n_times
y23 = epochs_MVP_R23.events[:, 2] # target
In [89]:
# define the Temporal Generalization object
time_gen = GeneralizingEstimator(clf, n_jobs=2, scoring=scoring_function1)
# 13
scores13 = cross_val_multiscore(time_gen, X13, y13, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scores13 = np.mean(scores13, axis=0)
# 12
scores12 = cross_val_multiscore(time_gen, X12, y12, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scores12 = np.mean(scores12, axis=0)
# 12
scores23 = cross_val_multiscore(time_gen, X23, y23, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scores23 = np.mean(scores23, axis=0)
In [96]:
fig, ax = plt.subplots(1,1)
im = ax.imshow(scores13, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_R13.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization R1 and R2')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
fig, ax = plt.subplots(1,1)
im = ax.imshow(scores12, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_R13.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization R1 and R2')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
fig, ax = plt.subplots(1,1)
im = ax.imshow(scores23, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_R23.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization R2 and R3')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
In [103]:
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
fig, ax = plt.subplots()
ax.plot(epochs_MVP_R13.times, np.diag(scores13), label='R1 vs R3')
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
ax.plot(epochs_MVP_R12.times, np.diag(scores12), label='R1 vs R2')
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
ax.plot(epochs_MVP_R23.times, np.diag(scores23), label='R2 vs R3')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC')
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Decoding MEG sensors over time Rings')
plt.show()
In [11]:
tmp_id_Q = np.arange(101,137,1).reshape(12, 3)
Q1_id = tmp_id_Q[[0,4,8],:].ravel()
Q2_id = tmp_id_Q[[1,5,9],:].ravel()
Q3_id = tmp_id_Q[[2,6,10],:].ravel()
Q4_id = tmp_id_Q[[3,7,11],:].ravel()
In [12]:
new_events_Q = mne.merge_events(events, Q1_id, 1)
new_events_Q = mne.merge_events(new_events_Q, Q2_id, 2)
new_events_Q = mne.merge_events(new_events_Q, Q3_id, 3)
new_events_Q = mne.merge_events(new_events_Q, Q4_id, 4)
In [15]:
picks = mne.pick_types(raw.info, meg=True, eog=False, misc=False)
new_event_Q12_id = dict(Q1=1, Q2=2)
epochs_MVP_Q12 = mne.Epochs(raw = raw, events = new_events_Q, event_id = new_event_Q12_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4, picks = picks)
new_event_Q13_id = dict(Q1=1, Q3=3)
epochs_MVP_Q13 = mne.Epochs(raw = raw, events = new_events_Q, event_id = new_event_Q13_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4, picks = picks)
new_event_Q14_id = dict(Q1=1, Q4=4)
epochs_MVP_Q14 = mne.Epochs(raw = raw, events = new_events_Q, event_id = new_event_Q14_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4, picks = picks)
new_event_Q23_id = dict(Q2=2, Q3=3)
epochs_MVP_Q23 = mne.Epochs(raw = raw, events = new_events_Q, event_id = new_event_Q23_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4, picks = picks)
new_event_Q24_id = dict(Q2=2, Q4=4)
epochs_MVP_Q24 = mne.Epochs(raw = raw, events = new_events_Q, event_id = new_event_Q24_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4, picks = picks)
new_event_Q34_id = dict(Q3=3, Q4=4)
epochs_MVP_Q34 = mne.Epochs(raw = raw, events = new_events_Q, event_id = new_event_Q34_id, tmin = tmin, tmax = tmax, proj=True,
baseline=None, detrend = 0, preload=True,
reject=reject, decim=4, picks = picks)
In [18]:
XQ12 = epochs_MVP_Q12.get_data() # MEG signals: n_epochs, n_channels, n_times
yQ12 = epochs_MVP_Q12.events[:, 2] # target
XQ13 = epochs_MVP_Q13.get_data() # MEG signals: n_epochs, n_channels, n_times
yQ13 = epochs_MVP_Q13.events[:, 2] # target
XQ14 = epochs_MVP_Q14.get_data() # MEG signals: n_epochs, n_channels, n_times
yQ14 = epochs_MVP_Q14.events[:, 2] # target
XQ23 = epochs_MVP_Q23.get_data() # MEG signals: n_epochs, n_channels, n_times
yQ23 = epochs_MVP_Q23.events[:, 2] # target
XQ24 = epochs_MVP_Q24.get_data() # MEG signals: n_epochs, n_channels, n_times
yQ24 = epochs_MVP_Q24.events[:, 2] # target
XQ34 = epochs_MVP_Q34.get_data() # MEG signals: n_epochs, n_channels, n_times
yQ34 = epochs_MVP_Q34.events[:, 2] # target
In [19]:
# define the Temporal Generalization object
clf = make_pipeline(StandardScaler(), LinearModel(LogisticRegression()))
time_gen = GeneralizingEstimator(clf, n_jobs=2, scoring=scoring_function1)
# 12
scoresQ12 = cross_val_multiscore(time_gen, XQ12, yQ12, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scoresQ12 = np.mean(scoresQ12, axis=0)
# 13
scoresQ13 = cross_val_multiscore(time_gen, XQ13, yQ13, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scoresQ13 = np.mean(scoresQ13, axis=0)
# 14
scoresQ14 = cross_val_multiscore(time_gen, XQ14, yQ14, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scoresQ14 = np.mean(scoresQ14, axis=0)
# 23
scoresQ23 = cross_val_multiscore(time_gen, XQ23, yQ23, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scoresQ23 = np.mean(scoresQ23, axis=0)
# 24
scoresQ24 = cross_val_multiscore(time_gen, XQ24, yQ24, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scoresQ24 = np.mean(scoresQ24, axis=0)
# 34
scoresQ34 = cross_val_multiscore(time_gen, XQ34, yQ34, cv=5, n_jobs=2)
# Mean scores across cross-validation splits
scoresQ34 = np.mean(scoresQ34, axis=0)
In [20]:
# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
fig, ax = plt.subplots()
ax.plot(epochs_MVP_Q12.times, np.diag(scoresQ12), label='Q1 vs Q2')
ax.plot(epochs_MVP_Q13.times, np.diag(scoresQ13), label='Q1 vs Q3')
ax.plot(epochs_MVP_Q14.times, np.diag(scoresQ14), label='Q1 vs Q4')
ax.plot(epochs_MVP_Q23.times, np.diag(scoresQ23), label='Q2 vs Q3')
ax.plot(epochs_MVP_Q24.times, np.diag(scoresQ24), label='Q2 vs Q4')
ax.plot(epochs_MVP_Q34.times, np.diag(scoresQ34), label='Q3 vs Q4')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC')
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Decoding MEG sensors over time in different quadrants')
plt.show()
In [21]:
fig, ax = plt.subplots(1,1)
im = ax.imshow(scoresQ14, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_Q14.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization Q1 and Q4')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
In [22]:
fig, ax = plt.subplots(1,1)
im = ax.imshow(scoresQ24, interpolation='lanczos', origin='lower', cmap='RdBu_r',
extent=epochs_MVP_Q24.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal Generalization Q2 and Q4')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)
plt.show()
In [ ]: